package com.pointi.emoEngine;

public class DSP_FFT
{
	DSP_FFT() {
	}

	public final double M_PI_2 = 1.5707963 ;
	public final int RDFT_LOOP_DIV = 64 ;
	public final int DCST_LOOP_DIV = 64 ;
	public void FFT(double [] sample, int sampleSize, int fft_size, double [] spectrum)
	{
		int i ;
		double [] K_temp_ = new double [(fft_size)*2] ; 

		//	Zero-Padding
		for (i=0 ; i<fft_size ; i++)
		{
			if (i < sampleSize)
			{
				K_temp_[i*2+1] = 0 ;	
				K_temp_[i*2] = sample[i] ;
			}
			else
			{
				K_temp_[i*2+1] = 0 ;	
				K_temp_[i*2] = 0 ;
			}
		}
		spectrum[0] = 0.0 ;
		cdft(fft_size*2,-1,K_temp_) ;

		// Compute magnitude 
		for (i=1 ; i < (fft_size / 2 + 1) ; i++)
		{
			spectrum[i] = Math.sqrt(K_temp_[2*i]*K_temp_[2*i] + K_temp_[2*i+1]*K_temp_[2*i+1]) ;
			spectrum[fft_size-i] = spectrum[i] ;
		}
		K_temp_ = null ;
	}

	public void divide(double [] x, int n)
	{
		int i ;
		double inv = (double)(1.0 / n) ;

		for(i=0 ; i < n ; i++)
		{
			x[2*i] = inv * x[2*i] ;
			x[2*i+1] = inv * x[2*i+1] ;
		}
	}

	public void cdft(int n, int isgn, double [] a)
	{
		if (n > 4)
		{
			if (isgn >= 0)
			{
				bitrv2(n, a) ;
				cftfsub(n, a) ;
			}
			else
			{
				bitrv2conj(n, a) ;
				cftbsub(n, a) ;
			}
		}
		else if (n == 4)
		{
			cftfsub(n, a) ;
		}
	}

	public void bitrv2(int n, double [] a)
	{
		int j0, k0, j1, k1, l, m, i, j, k ;
		double xr, xi, yr, yi ;

		l = n >> 2 ;
			m = 2 ;
			while (m < l)
			{
				l >>= 1 ;
			m <<= 1 ;
			}
			if (m == l) {
				j0 = 0 ;
				for (k0 = 0 ; k0 < m ; k0 += 2) {
					k = k0 ;
					for (j = j0 ; j < j0 + k0 ; j += 2) {
						xr = a[j] ;
						xi = a[j + 1] ;
						yr = a[k] ;
						yi = a[k + 1] ;
						a[j] = yr ;
						a[j + 1] = yi ;
						a[k] = xr ;
						a[k + 1] = xi ;
						j1 = j + m ;
						k1 = k + 2 * m ;
						xr = a[j1] ;
						xi = a[j1 + 1] ;
						yr = a[k1] ;
						yi = a[k1 + 1] ;
						a[j1] = yr ;
						a[j1 + 1] = yi ;
						a[k1] = xr ;
						a[k1 + 1] = xi ;
						j1 += m ;
						k1 -= m ;
						xr = a[j1] ;
						xi = a[j1 + 1] ;
						yr = a[k1] ;
						yi = a[k1 + 1] ;
						a[j1] = yr ;
						a[j1 + 1] = yi ;
						a[k1] = xr ;
						a[k1 + 1] = xi ;
						j1 += m ;
						k1 += 2 * m ;
						xr = a[j1] ;
						xi = a[j1 + 1] ;
						yr = a[k1] ;
						yi = a[k1 + 1] ;
						a[j1] = yr ;
						a[j1 + 1] = yi ;
						a[k1] = xr ;
						a[k1 + 1] = xi ;
						for (i = n >> 1 ; i > (k ^= i) ; i >>= 1) ;
					}
					j1 = j0 + k0 + m ;
					k1 = j1 + m ;
					xr = a[j1] ;
					xi = a[j1 + 1] ;
					yr = a[k1] ;
					yi = a[k1 + 1] ;
					a[j1] = yr ;
					a[j1 + 1] = yi ;
					a[k1] = xr ;
					a[k1 + 1] = xi ;
					for (i = n >> 1 ; i > (j0 ^= i) ; i >>= 1) ;
				}
			} else {
				j0 = 0 ;
				for (k0 = 2 ; k0 < m ; k0 += 2) {
					for (i = n >> 1 ; i > (j0 ^= i) ; i >>= 1) ;
					k = k0 ;
					for (j = j0 ; j < j0 + k0 ; j += 2) {
						xr = a[j] ;
						xi = a[j + 1] ;
						yr = a[k] ;
						yi = a[k + 1] ;
						a[j] = yr ;
						a[j + 1] = yi ;
						a[k] = xr ;
						a[k + 1] = xi ;
						j1 = j + m ;
						k1 = k + m ;
						xr = a[j1] ;
						xi = a[j1 + 1] ;
						yr = a[k1] ;
						yi = a[k1 + 1] ;
						a[j1] = yr ;
						a[j1 + 1] = yi ;
						a[k1] = xr ;
						a[k1 + 1] = xi ;
						for (i = n >> 1 ; i > (k ^= i) ; i >>= 1) ;
					}
				}
			}
	}

	public void bitrv2conj(int n, double [] a)
	{
		int j0, k0, j1, k1, l, m, i, j, k ;
		double xr, xi, yr, yi ;

		l = n >> 2 ;
		m = 2 ;
		while (m < l) {
			l >>= 1 ;
		m <<= 1 ;
		}
		if (m == l) {
			j0 = 0 ;
			for (k0 = 0 ; k0 < m ; k0 += 2) {
				k = k0 ;
				for (j = j0 ; j < j0 + k0 ; j += 2) {
					xr = a[j] ;
					xi = -a[j + 1] ;
					yr = a[k] ;
					yi = -a[k + 1] ;
					a[j] = yr ;
					a[j + 1] = yi ;
					a[k] = xr ;
					a[k + 1] = xi ;
					j1 = j + m ;
					k1 = k + 2 * m ;
					xr = a[j1] ;
					xi = -a[j1 + 1] ;
					yr = a[k1] ;
					yi = -a[k1 + 1] ;
					a[j1] = yr ;
					a[j1 + 1] = yi ;
					a[k1] = xr ;
					a[k1 + 1] = xi ;
					j1 += m ;
					k1 -= m ;
					xr = a[j1] ;
					xi = -a[j1 + 1] ;
					yr = a[k1] ;
					yi = -a[k1 + 1] ;
					a[j1] = yr ;
					a[j1 + 1] = yi ;
					a[k1] = xr ;
					a[k1 + 1] = xi ;
					j1 += m ;
					k1 += 2 * m ;
					xr = a[j1] ;
					xi = -a[j1 + 1] ;
					yr = a[k1] ;
					yi = -a[k1 + 1] ;
					a[j1] = yr ;
					a[j1 + 1] = yi ;
					a[k1] = xr ;
					a[k1 + 1] = xi ;
					for (i = n >> 1 ; i > (k ^= i) ; i >>= 1) ;
				}
				k1 = j0 + k0 ;
				a[k1 + 1] = -a[k1 + 1] ;
				j1 = k1 + m ;
				k1 = j1 + m ;
				xr = a[j1] ;
				xi = -a[j1 + 1] ;
				yr = a[k1] ;
				yi = -a[k1 + 1] ;
				a[j1] = yr ;
				a[j1 + 1] = yi ;
				a[k1] = xr ;
				a[k1 + 1] = xi ;
				k1 += m ;
				a[k1 + 1] = -a[k1 + 1] ;
				for (i = n >> 1 ; i > (j0 ^= i) ; i >>= 1) ;
			}
		} else {
			a[1] = -a[1] ;
			a[m + 1] = -a[m + 1] ;
			j0 = 0 ;
			for (k0 = 2 ; k0 < m ; k0 += 2) {
				for (i = n >> 1 ; i > (j0 ^= i) ; i >>= 1) ;
				k = k0 ;
				for (j = j0 ; j < j0 + k0 ; j += 2) {
					xr = a[j] ;
					xi = -a[j + 1] ;
					yr = a[k] ;
					yi = -a[k + 1] ;
					a[j] = yr ;
					a[j + 1] = yi ;
					a[k] = xr ;
					a[k + 1] = xi ;
					j1 = j + m ;
					k1 = k + m ;
					xr = a[j1] ;
					xi = -a[j1 + 1] ;
					yr = a[k1] ;
					yi = -a[k1 + 1] ;
					a[j1] = yr ;
					a[j1 + 1] = yi ;
					a[k1] = xr ;
					a[k1 + 1] = xi ;
					for (i = n >> 1 ; i > (k ^= i) ; i >>= 1) ;
				}
				k1 = j0 + k0 ;
				a[k1 + 1] = -a[k1 + 1] ;
				a[k1 + m + 1] = -a[k1 + m + 1] ;
			}
		}
	}

	public void bitrv1(int n, double [] a)
	{
		int j0, k0, j1, k1, l, m, i, j, k ;
		double x ;

		l = n >> 2 ;
		m = 1 ;
		while (m < l)
		{
			l >>= 1 ;
			m <<= 1 ;
		}
		if (m == l) {
			j0 = 0 ;
			for (k0 = 0 ; k0 < m ; k0++)
			{
				k = k0 ;
				for (j = j0 ; j < j0 + k0 ; j++)
				{
					x = a[j] ;
					a[j] = a[k] ;
					a[k] = x ;
					j1 = j + m ;
					k1 = k + 2 * m ;
					x = a[j1] ;
					a[j1] = a[k1] ;
					a[k1] = x ;
					j1 += m ;
					k1 -= m ;
					x = a[j1] ;
					a[j1] = a[k1] ;
					a[k1] = x ;
					j1 += m ;
					k1 += 2 * m ;
					x = a[j1] ;
					a[j1] = a[k1] ;
					a[k1] = x ;
					for (i = n >> 1 ; i > (k ^= i) ; i >>= 1) ;
				}
				j1 = j0 + k0 + m ;
				k1 = j1 + m ;
				x = a[j1] ;
				a[j1] = a[k1] ;
				a[k1] = x ;
				for (i = n >> 1 ; i > (j0 ^= i) ; i >>= 1) ;
			}
		}
		else
		{
			j0 = 0 ;
			for (k0 = 1 ; k0 < m ; k0++)
			{
				for (i = n >> 1 ; i > (j0 ^= i) ; i >>= 1) ;
				k = k0 ;
				for (j = j0 ; j < j0 + k0 ; j++)
				{
					x = a[j] ;
					a[j] = a[k] ;
					a[k] = x ;
					j1 = j + m ;
					k1 = k + m ;
					x = a[j1] ;
					a[j1] = a[k1] ;
					a[k1] = x ;
					for (i = n >> 1 ; i > (k ^= i) ; i >>= 1) ;
				}
			}
		}
	}

	public void cftfsub(int n, double [] a)
	{
		int j, j1, j2, j3, l ;
		double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i ;

		l = 2 ;
		if (n > 8) {
			cft1st(n, a) ;
			l = 8 ;
			while ((l << 2) < n) {
				cftmdl(n, l, a) ;
				l <<= 2 ;
			}
		}
		if ((l << 2) == n) {
			for (j = 0 ; j < l ; j += 2) {
				j1 = j + l ;
				j2 = j1 + l ;
				j3 = j2 + l ;
				x0r = a[j] + a[j1] ;
				x0i = a[j + 1] + a[j1 + 1] ;
				x1r = a[j] - a[j1] ;
				x1i = a[j + 1] - a[j1 + 1] ;
				x2r = a[j2] + a[j3] ;
				x2i = a[j2 + 1] + a[j3 + 1] ;
				x3r = a[j2] - a[j3] ;
				x3i = a[j2 + 1] - a[j3 + 1] ;
				a[j] = x0r + x2r ;
				a[j + 1] = x0i + x2i ;
				a[j2] = x0r - x2r ;
				a[j2 + 1] = x0i - x2i ;
				a[j1] = x1r - x3i ;
				a[j1 + 1] = x1i + x3r ;
				a[j3] = x1r + x3i ;
				a[j3 + 1] = x1i - x3r ;
			}
		} else {
			for (j = 0 ; j < l ; j += 2) {
				j1 = j + l ;
				x0r = a[j] - a[j1] ;
				x0i = a[j + 1] - a[j1 + 1] ;
				a[j] += a[j1] ;
				a[j + 1] += a[j1 + 1] ;
				a[j1] = x0r ;
				a[j1 + 1] = x0i ;
			}
		}
	}

	public void cftbsub(int n, double [] a)
	{
		int j, j1, j2, j3, l ;
		double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i ;

		l = 2 ;
		if (n > 8) {
			cft1st(n, a) ;
			l = 8 ;
			while ((l << 2) < n) {
				cftmdl(n, l, a) ;
				l <<= 2 ;
			}
		}
		if ((l << 2) == n) {
			for (j = 0 ; j < l ; j += 2) {
				j1 = j + l ;
				j2 = j1 + l ;
				j3 = j2 + l ;
				x0r = a[j] + a[j1] ;
				x0i = -a[j + 1] - a[j1 + 1] ;
				x1r = a[j] - a[j1] ;
				x1i = -a[j + 1] + a[j1 + 1] ;
				x2r = a[j2] + a[j3] ;
				x2i = a[j2 + 1] + a[j3 + 1] ;
				x3r = a[j2] - a[j3] ;
				x3i = a[j2 + 1] - a[j3 + 1] ;
				a[j] = x0r + x2r ;
				a[j + 1] = x0i - x2i ;
				a[j2] = x0r - x2r ;
				a[j2 + 1] = x0i + x2i ;
				a[j1] = x1r - x3i ;
				a[j1 + 1] = x1i - x3r ;
				a[j3] = x1r + x3i ;
				a[j3 + 1] = x1i + x3r ;
			}
		} else {
			for (j = 0 ; j < l ; j += 2) {
				j1 = j + l ;
				x0r = a[j] - a[j1] ;
				x0i = -a[j + 1] + a[j1 + 1] ;
				a[j] += a[j1] ;
				a[j + 1] = -a[j + 1] - a[j1 + 1] ;
				a[j1] = x0r ;
				a[j1 + 1] = x0i ;
			}
		}
	}

	public void cft1st(int n, double [] a)
	{
		int j, kj, kr ;
		double ew, wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i ;
		double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i ;

		x0r = a[0] + a[2] ;
		x0i = a[1] + a[3] ;
		x1r = a[0] - a[2] ;
		x1i = a[1] - a[3] ;
		x2r = a[4] + a[6] ;
		x2i = a[5] + a[7] ;
		x3r = a[4] - a[6] ;
		x3i = a[5] - a[7] ;
		a[0] = x0r + x2r ;
		a[1] = x0i + x2i ;
		a[4] = x0r - x2r ;
		a[5] = x0i - x2i ;
		a[2] = x1r - x3i ;
		a[3] = x1i + x3r ;
		a[6] = x1r + x3i ;
		a[7] = x1i - x3r ;
		wn4r = (double)Math.cos(M_PI_2 * 0.5) ;
		x0r = a[8] + a[10] ;
		x0i = a[9] + a[11] ;
		x1r = a[8] - a[10] ;
		x1i = a[9] - a[11] ;
		x2r = a[12] + a[14] ;
		x2i = a[13] + a[15] ;
		x3r = a[12] - a[14] ;
		x3i = a[13] - a[15] ;
		a[8] = x0r + x2r ;
		a[9] = x0i + x2i ;
		a[12] = x2i - x0i ;
		a[13] = x0r - x2r ;
		x0r = x1r - x3i ;
		x0i = x1i + x3r ;
		a[10] = wn4r * (x0r - x0i) ;
		a[11] = wn4r * (x0r + x0i) ;
		x0r = x3i + x1r ;
		x0i = x3r - x1i ;
		a[14] = wn4r * (x0i - x0r) ;
		a[15] = wn4r * (x0i + x0r) ;
		ew = M_PI_2 / n ;
		kr = 0 ;
		for (j = 16 ; j < n ; j += 16) {
			for (kj = n >> 2 ; kj > (kr ^= kj) ; kj >>= 1) ;
			wk1r = (double)Math.cos(ew * kr) ;
			wk1i = (double)Math.sin(ew * kr) ;
			wk2r = 1 - 2 * wk1i * wk1i ;
			wk2i = 2 * wk1i * wk1r ;
			wk3r = wk1r - 2 * wk2i * wk1i ;
			wk3i = 2 * wk2i * wk1r - wk1i ;
			x0r = a[j] + a[j + 2] ;
			x0i = a[j + 1] + a[j + 3] ;
			x1r = a[j] - a[j + 2] ;
			x1i = a[j + 1] - a[j + 3] ;
			x2r = a[j + 4] + a[j + 6] ;
			x2i = a[j + 5] + a[j + 7] ;
			x3r = a[j + 4] - a[j + 6] ;
			x3i = a[j + 5] - a[j + 7] ;
			a[j] = x0r + x2r ;
			a[j + 1] = x0i + x2i ;
			x0r -= x2r ;
			x0i -= x2i ;
			a[j + 4] = wk2r * x0r - wk2i * x0i ;
			a[j + 5] = wk2r * x0i + wk2i * x0r ;
			x0r = x1r - x3i ;
			x0i = x1i + x3r ;
			a[j + 2] = wk1r * x0r - wk1i * x0i ;
			a[j + 3] = wk1r * x0i + wk1i * x0r ;
			x0r = x1r + x3i ;
			x0i = x1i - x3r ;
			a[j + 6] = wk3r * x0r - wk3i * x0i ;
			a[j + 7] = wk3r * x0i + wk3i * x0r ;
			x0r = wn4r * (wk1r - wk1i) ;
			wk1i = wn4r * (wk1r + wk1i) ;
			wk1r = x0r ;
			wk3r = wk1r - 2 * wk2r * wk1i ;
			wk3i = 2 * wk2r * wk1r - wk1i ;
			x0r = a[j + 8] + a[j + 10] ;
			x0i = a[j + 9] + a[j + 11] ;
			x1r = a[j + 8] - a[j + 10] ;
			x1i = a[j + 9] - a[j + 11] ;
			x2r = a[j + 12] + a[j + 14] ;
			x2i = a[j + 13] + a[j + 15] ;
			x3r = a[j + 12] - a[j + 14] ;
			x3i = a[j + 13] - a[j + 15] ;
			a[j + 8] = x0r + x2r ;
			a[j + 9] = x0i + x2i ;
			x0r -= x2r ;
			x0i -= x2i ;
			a[j + 12] = -wk2i * x0r - wk2r * x0i ;
			a[j + 13] = -wk2i * x0i + wk2r * x0r ;
			x0r = x1r - x3i ;
			x0i = x1i + x3r ;
			a[j + 10] = wk1r * x0r - wk1i * x0i ;
			a[j + 11] = wk1r * x0i + wk1i * x0r ;
			x0r = x1r + x3i ;
			x0i = x1i - x3r ;
			a[j + 14] = wk3r * x0r - wk3i * x0i ;
			a[j + 15] = wk3r * x0i + wk3i * x0r ;
		}
	}

	public void cftmdl(int n, int l, double [] a)
	{
		int j, j1, j2, j3, k, kj, kr, m, m2 ;
		double ew, wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i ;
		double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i ;

		m = l << 2 ;
		for (j = 0 ; j < l ; j += 2) {
			j1 = j + l ;
			j2 = j1 + l ;
			j3 = j2 + l ;
			x0r = a[j] + a[j1] ;
			x0i = a[j + 1] + a[j1 + 1] ;
			x1r = a[j] - a[j1] ;
			x1i = a[j + 1] - a[j1 + 1] ;
			x2r = a[j2] + a[j3] ;
			x2i = a[j2 + 1] + a[j3 + 1] ;
			x3r = a[j2] - a[j3] ;
			x3i = a[j2 + 1] - a[j3 + 1] ;
			a[j] = x0r + x2r ;
			a[j + 1] = x0i + x2i ;
			a[j2] = x0r - x2r ;
			a[j2 + 1] = x0i - x2i ;
			a[j1] = x1r - x3i ;
			a[j1 + 1] = x1i + x3r ;
			a[j3] = x1r + x3i ;
			a[j3 + 1] = x1i - x3r ;
		}
		wn4r = (double)Math.cos(M_PI_2 * 0.5) ;
		for (j = m ; j < l + m ; j += 2) {
			j1 = j + l ;
			j2 = j1 + l ;
			j3 = j2 + l ;
			x0r = a[j] + a[j1] ;
			x0i = a[j + 1] + a[j1 + 1] ;
			x1r = a[j] - a[j1] ;
			x1i = a[j + 1] - a[j1 + 1] ;
			x2r = a[j2] + a[j3] ;
			x2i = a[j2 + 1] + a[j3 + 1] ;
			x3r = a[j2] - a[j3] ;
			x3i = a[j2 + 1] - a[j3 + 1] ;
			a[j] = x0r + x2r ;
			a[j + 1] = x0i + x2i ;
			a[j2] = x2i - x0i ;
			a[j2 + 1] = x0r - x2r ;
			x0r = x1r - x3i ;
			x0i = x1i + x3r ;
			a[j1] = wn4r * (x0r - x0i) ;
			a[j1 + 1] = wn4r * (x0r + x0i) ;
			x0r = x3i + x1r ;
			x0i = x3r - x1i ;
			a[j3] = wn4r * (x0i - x0r) ;
			a[j3 + 1] = wn4r * (x0i + x0r) ;
		}
		ew = M_PI_2 / n ;
		kr = 0 ;
		m2 = 2 * m ;
		for (k = m2 ; k < n ; k += m2) {
			for (kj = n >> 2 ; kj > (kr ^= kj) ; kj >>= 1) ;
			wk1r = (double)Math.cos(ew * kr) ;
			wk1i = (double)Math.sin(ew * kr) ;
			wk2r = 1 - 2 * wk1i * wk1i ;
			wk2i = 2 * wk1i * wk1r ;
			wk3r = wk1r - 2 * wk2i * wk1i ;
			wk3i = 2 * wk2i * wk1r - wk1i ;
			for (j = k ; j < l + k ; j += 2) {
				j1 = j + l ;
				j2 = j1 + l ;
				j3 = j2 + l ;
				x0r = a[j] + a[j1] ;
				x0i = a[j + 1] + a[j1 + 1] ;
				x1r = a[j] - a[j1] ;
				x1i = a[j + 1] - a[j1 + 1] ;
				x2r = a[j2] + a[j3] ;
				x2i = a[j2 + 1] + a[j3 + 1] ;
				x3r = a[j2] - a[j3] ;
				x3i = a[j2 + 1] - a[j3 + 1] ;
				a[j] = x0r + x2r ;
				a[j + 1] = x0i + x2i ;
				x0r -= x2r ;
				x0i -= x2i ;
				a[j2] = wk2r * x0r - wk2i * x0i ;
				a[j2 + 1] = wk2r * x0i + wk2i * x0r ;
				x0r = x1r - x3i ;
				x0i = x1i + x3r ;
				a[j1] = wk1r * x0r - wk1i * x0i ;
				a[j1 + 1] = wk1r * x0i + wk1i * x0r ;
				x0r = x1r + x3i ;
				x0i = x1i - x3r ;
				a[j3] = wk3r * x0r - wk3i * x0i ;
				a[j3 + 1] = wk3r * x0i + wk3i * x0r ;
			}
			x0r = wn4r * (wk1r - wk1i) ;
			wk1i = wn4r * (wk1r + wk1i) ;
			wk1r = x0r ;
			wk3r = wk1r - 2 * wk2r * wk1i ;
			wk3i = 2 * wk2r * wk1r - wk1i ;
			for (j = k + m ; j < l + (k + m) ; j += 2) {
				j1 = j + l ;
				j2 = j1 + l ;
				j3 = j2 + l ;
				x0r = a[j] + a[j1] ;
				x0i = a[j + 1] + a[j1 + 1] ;
				x1r = a[j] - a[j1] ;
				x1i = a[j + 1] - a[j1 + 1] ;
				x2r = a[j2] + a[j3] ;
				x2i = a[j2 + 1] + a[j3 + 1] ;
				x3r = a[j2] - a[j3] ;
				x3i = a[j2 + 1] - a[j3 + 1] ;
				a[j] = x0r + x2r ;
				a[j + 1] = x0i + x2i ;
				x0r -= x2r ;
				x0i -= x2i ;
				a[j2] = -wk2i * x0r - wk2r * x0i ;
				a[j2 + 1] = -wk2i * x0i + wk2r * x0r ;
				x0r = x1r - x3i ;
				x0i = x1i + x3r ;
				a[j1] = wk1r * x0r - wk1i * x0i ;
				a[j1 + 1] = wk1r * x0i + wk1i * x0r ;
				x0r = x1r + x3i ;
				x0i = x1i - x3r ;
				a[j3] = wk3r * x0r - wk3i * x0i ;
				a[j3 + 1] = wk3r * x0i + wk3i * x0r ;
			}
		}
	}

	public void rftfsub(int n, double [] a)
	{
		int i, j, k ;
		double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi ;

		ec = 2 * M_PI_2 / n ;
		wkr = 0 ;
		wki = 0 ;
		wdi = (double)Math.cos(ec) ;
		wdr = (double)Math.sin(ec) ;
		wdi *= wdr ;
		wdr *= wdr ;
		w1r = 1 - 2 * wdr ;
		w1i = 2 * wdi ;
		ss = 2 * w1i ;
		j = n >> 1 ;
			while (j > 4 * RDFT_LOOP_DIV) {
				for (i = 0 ; i < RDFT_LOOP_DIV ; i++) {
					j -= 4 ;
					k = n - j ;
					xr = a[j + 2] - a[k - 2] ;
					xi = a[j + 3] + a[k - 1] ;
					yr = wdr * xr - wdi * xi ;
					yi = wdr * xi + wdi * xr ;
					a[j + 2] -= yr ;
					a[j + 3] -= yi ;
					a[k - 2] += yr ;
					a[k - 1] -= yi ;
					wkr += ss * wdi ;
					wki += ss * (0.5 - wdr) ;
					xr = a[j] - a[k] ;
					xi = a[j + 1] + a[k + 1] ;
					yr = wkr * xr - wki * xi ;
					yi = wkr * xi + wki * xr ;
					a[j] -= yr ;
					a[j + 1] -= yi ;
					a[k] += yr ;
					a[k + 1] -= yi ;
					wdr += ss * wki ;
					wdi += ss * (0.5 - wkr) ;
				}
				wkr = 0.5 * (double)Math.sin(ec * j) ;
				wki = 0.5 * (double)Math.cos(ec * j) ;
				wdr = 0.5 - (wkr * w1r - wki * w1i) ;
				wdi = wkr * w1i + wki * w1r ;
				wkr = 0.5 - wkr ;
			}
			while (j > 4) {
				j -= 4 ;
				k = n - j ;
				xr = a[j + 2] - a[k - 2] ;
				xi = a[j + 3] + a[k - 1] ;
				yr = wdr * xr - wdi * xi ;
				yi = wdr * xi + wdi * xr ;
				a[j + 2] -= yr ;
				a[j + 3] -= yi ;
				a[k - 2] += yr ;
				a[k - 1] -= yi ;
				wkr += ss * wdi ;
				wki += ss * (0.5 - wdr) ;
				xr = a[j] - a[k] ;
				xi = a[j + 1] + a[k + 1] ;
				yr = wkr * xr - wki * xi ;
				yi = wkr * xi + wki * xr ;
				a[j] -= yr ;
				a[j + 1] -= yi ;
				a[k] += yr ;
				a[k + 1] -= yi ;
				wdr += ss * wki ;
				wdi += ss * (0.5 - wkr) ;
			}
			xr = a[2] - a[n - 2] ;
			xi = a[3] + a[n - 1] ;
			yr = wdr * xr - wdi * xi ;
			yi = wdr * xi + wdi * xr ;
			a[2] -= yr ;
			a[3] -= yi ;
			a[n - 2] += yr ;
			a[n - 1] -= yi ;
	}

	public void rftbsub(int n, double [] a)
	{
		int i, j, k ;
		double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi ;

		ec = 2 * M_PI_2 / n ;
		wkr = 0 ;
		wki = 0 ;
		wdi = (double)Math.cos(ec) ;
		wdr = (double)Math.sin(ec) ;
		wdi *= wdr ;
		wdr *= wdr ;
		w1r = 1 - 2 * wdr ;
		w1i = 2 * wdi ;
		ss = 2 * w1i ;
		j = n >> 1 ;
				a[j + 1] = -a[j + 1] ;
				while (j > 4 * RDFT_LOOP_DIV) {
					for (i = 0 ; i < RDFT_LOOP_DIV ; i++) {
						j -= 4 ;
						k = n - j ;
						xr = a[j + 2] - a[k - 2] ;
						xi = a[j + 3] + a[k - 1] ;
						yr = wdr * xr + wdi * xi ;
						yi = wdr * xi - wdi * xr ;
						a[j + 2] -= yr ;
						a[j + 3] = yi - a[j + 3] ;
						a[k - 2] += yr ;
						a[k - 1] = yi - a[k - 1] ;
						wkr += ss * wdi ;
						wki += ss * (0.5 - wdr) ;
						xr = a[j] - a[k] ;
						xi = a[j + 1] + a[k + 1] ;
						yr = wkr * xr + wki * xi ;
						yi = wkr * xi - wki * xr ;
						a[j] -= yr ;
						a[j + 1] = yi - a[j + 1] ;
						a[k] += yr ;
						a[k + 1] = yi - a[k + 1] ;
						wdr += ss * wki ;
						wdi += ss * (0.5 - wkr) ;
					}
					wkr = 0.5 * (double)Math.sin(ec * j) ;
					wki = 0.5 * (double)Math.cos(ec * j) ;
					wdr = 0.5 - (wkr * w1r - wki * w1i) ;
					wdi = wkr * w1i + wki * w1r ;
					wkr = 0.5 - wkr ;
				}
				while (j > 4) {
					j -= 4 ;
					k = n - j ;
					xr = a[j + 2] - a[k - 2] ;
					xi = a[j + 3] + a[k - 1] ;
					yr = wdr * xr + wdi * xi ;
					yi = wdr * xi - wdi * xr ;
					a[j + 2] -= yr ;
					a[j + 3] = yi - a[j + 3] ;
					a[k - 2] += yr ;
					a[k - 1] = yi - a[k - 1] ;
					wkr += ss * wdi ;
					wki += ss * (0.5 - wdr) ;
					xr = a[j] - a[k] ;
					xi = a[j + 1] + a[k + 1] ;
					yr = wkr * xr + wki * xi ;
					yi = wkr * xi - wki * xr ;
					a[j] -= yr ;
					a[j + 1] = yi - a[j + 1] ;
					a[k] += yr ;
					a[k + 1] = yi - a[k + 1] ;
					wdr += ss * wki ;
					wdi += ss * (0.5 - wkr) ;
				}
				xr = a[2] - a[n - 2] ;
				xi = a[3] + a[n - 1] ;
				yr = wdr * xr + wdi * xi ;
				yi = wdr * xi - wdi * xr ;
				a[2] -= yr ;
				a[3] = yi - a[3] ;
				a[n - 2] += yr ;
				a[n - 1] = yi - a[n - 1] ;
				a[1] = -a[1] ;
	}

	public void dctsub(int n, double [] a)
	{
		int i, j, k, m ;
		double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi ;

		ec = M_PI_2 / n ;
		wkr = 0.5 ;
		wki = 0.5 ;
		w1r = (double)Math.cos(ec) ;
		w1i = (double)Math.sin(ec) ;
		wdr = 0.5 * (w1r - w1i) ;
		wdi = 0.5 * (w1r + w1i) ;
		ss = 2 * w1i ;
		m = n >> 1 ;
			j = 0 ;
			while (j < m - 2 * DCST_LOOP_DIV) {
				for (i = 0 ; i < DCST_LOOP_DIV ; i++) {
					j += 2 ;
					k = n - j ;
					xr = wdi * a[j - 1] - wdr * a[k + 1] ;
					xi = wdr * a[j - 1] + wdi * a[k + 1] ;
					wkr -= ss * wdi ;
					wki += ss * wdr ;
					yr = wki * a[j] - wkr * a[k] ;
					yi = wkr * a[j] + wki * a[k] ;
					wdr -= ss * wki ;
					wdi += ss * wkr ;
					a[k + 1] = xr ;
					a[k] = yr ;
					a[j - 1] = xi ;
					a[j] = yi ;
				}
				wdr = (double)Math.cos(ec * j) ;
				wdi = (double)Math.sin(ec * j) ;
				wkr = 0.5 * (wdr - wdi) ;
				wki = 0.5 * (wdr + wdi) ;
				wdr = wkr * w1r - wki * w1i ;
				wdi = wkr * w1i + wki * w1r ;
			}
			while (j < m - 2) {
				j += 2 ;
				k = n - j ;
				xr = wdi * a[j - 1] - wdr * a[k + 1] ;
				xi = wdr * a[j - 1] + wdi * a[k + 1] ;
				wkr -= ss * wdi ;
				wki += ss * wdr ;
				yr = wki * a[j] - wkr * a[k] ;
				yi = wkr * a[j] + wki * a[k] ;
				wdr -= ss * wki ;
				wdi += ss * wkr ;
				a[k + 1] = xr ;
				a[k] = yr ;
				a[j - 1] = xi ;
				a[j] = yi ;
			}
			xr = wdi * a[m - 1] - wdr * a[m + 1] ;
			a[m - 1] = wdr * a[m - 1] + wdi * a[m + 1] ;
			a[m + 1] = xr ;
			a[m] *= wki + ss * wdr ;
	}

	public void dstsub(int n, double [] a)
	{
		int i, j, k, m ;
		double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi ;
		ec = M_PI_2 / n ;
		wkr = 0.5 ;
		wki = 0.5 ;
		w1r = (double)Math.cos(ec) ;
		w1i = (double)Math.sin(ec) ;
		wdr = 0.5 * (w1r - w1i) ;
		wdi = 0.5 * (w1r + w1i) ;
		ss = 2 * w1i ;
		m = n >> 1 ;
			j = 0 ;
			while (j < m - 2 * DCST_LOOP_DIV) {
				for (i = 0 ; i < DCST_LOOP_DIV ; i++) {
					j += 2 ;
					k = n - j ;
					xr = wdi * a[k + 1] - wdr * a[j - 1] ;
					xi = wdr * a[k + 1] + wdi * a[j - 1] ;
					wkr -= ss * wdi ;
					wki += ss * wdr ;
					yr = wki * a[k] - wkr * a[j] ;
					yi = wkr * a[k] + wki * a[j] ;
					wdr -= ss * wki ;
					wdi += ss * wkr ;
					a[j - 1] = xr ;
					a[j] = yr ;
					a[k + 1] = xi ;
					a[k] = yi ;
				}
				wdr = (double)Math.cos(ec * j) ;
				wdi = (double)Math.sin(ec * j) ;
				wkr = 0.5 * (wdr - wdi) ;
				wki = 0.5 * (wdr + wdi) ;
				wdr = wkr * w1r - wki * w1i ;
				wdi = wkr * w1i + wki * w1r ;
			}
			while (j < m - 2) {
				j += 2 ;
				k = n - j ;
				xr = wdi * a[k + 1] - wdr * a[j - 1] ;
				xi = wdr * a[k + 1] + wdi * a[j - 1] ;
				wkr -= ss * wdi ;
				wki += ss * wdr ;
				yr = wki * a[k] - wkr * a[j] ;
				yi = wkr * a[k] + wki * a[j] ;
				wdr -= ss * wki ;
				wdi += ss * wkr ;
				a[j - 1] = xr ;
				a[j] = yr ;
				a[k + 1] = xi ;
				a[k] = yi ;
			}
			xr = wdi * a[m + 1] - wdr * a[m - 1] ;
			a[m + 1] = wdr * a[m + 1] + wdi * a[m - 1] ;
			a[m - 1] = xr ;
			a[m] *= wki + ss * wdr ;
	}

	public void dctsub4(int n, double [] a)
	{
		int m ;
		double wki, wdr, wdi, xr ;

		wki = (double)Math.cos(M_PI_2 * 0.5) ;
		m = n >> 1 ;
			if (m == 2) {
				wdr = wki * Math.sin(M_PI_2 * 0.25) ;
				wdi = wki * Math.cos(M_PI_2 * 0.25) ;
				xr = wdi * a[1] - wdr * a[3] ;
				a[1] = wdr * a[1] + wdi * a[3] ;
				a[3] = xr ;
			}
			a[m] *= wki ;
	}

	public void dstsub4(int n, double [] a)
	{
		int m ;
		double wki, wdr, wdi, xr ;

		wki = Math.cos(M_PI_2 * 0.5) ;
		m = n >> 1 ;
			if (m == 2) {
				wdr = wki * Math.sin(M_PI_2 * 0.25) ;
				wdi = wki * Math.cos(M_PI_2 * 0.25) ;
				xr = wdi * a[3] - wdr * a[1] ;
				a[3] = wdr * a[3] + wdi * a[1] ;
				a[1] = xr ;
			}
			a[m] *= wki ;
	}
}